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A method for computing discrete-time state-space models of linearized unsteady 
aerodynamic behavior directly from aeroelastic CFD codes is presented. The method 
involves the treatment of CFD-based pulse responses as Markov parameters for use in 
a system identification/realization algorithm. Results are presented for the AGARD 
445.6 Aeroelastic Wing with four aeroelastic modes at a Mach number of 0.96 us- 
ing the EZNSS Euler/Navier-Stokes flow solver with aeroelastic capability. The Sys- 
tem/Observer/Controller Identification Toolbox (SOCIT) algorithm, based on the Ho- 
Kalman realization algorithm, is used to generate 15th- and 32nd-order discrete-time 
state-space models of the unsteady aerodynamic response of the wing over the entire 
frequency range of interest. 


Introduction 

T HE inclusion of CFD-based analyses into disci- 
plines such as aeroelasticity, aeroservoelasticity, 
and optimization is currently not performed on a rou- 
tine basis due to the high computational costs of the 
CFD portion of the analyses. One solution to this 
problem is the development of CFD-based reduced- 
order models (ROMs). These ROMs capture the 
essence of the dynamical system under investigation 
while reducing the complexity of the computational 
model. 

At present, the development of CFD-based ROMs 
is an area of active research at several industry, gov- 
ernment, and academic institutions . 1 Development of 
ROMs based on the Volterra theory is one of sev- 
eral ROM methods currently under development . 2 " 5 
Reduced-order models based on the Volterra theory 
have been applied successfully to Euler and Navier- 
Stokes models of nonlinear unsteady aerodynamic and 
aeroelastic systems. Volterra-based ROMs are based 
on the creation of unsteady aerodynamic impulse or 
step responses that are then used in a convolution 
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scheme to provide the linearized and nonlinear re- 
sponses of the system to arbitrary inputs. In this 
setting, the impulse or step responses ore the func- 
tional ROMs. 

A traditional approach for obtaining linearized gen- 
eralized aerodynamic forces (GAFs) from an aeroelas- 
tic CFD model is to perform a time-domain perturba- 
tion of all the modes (one mode at a time) resulting in 
a GAF influence coefficient matrix. This time-domain 
influence coefficient matrix is then transformed into 
the frequency domain 6 using standard Fourier trans- 
form techniques. The resultant frequency-domain 
GAFs can then be used in standard frequency- 
domain aeroelastic analyses. In addition, if time- 
domain aeroservoelastic (ASE) analyses are desired, 
the frequency-domain GAFs are transformed back into 
the time domain using rational function approxima- 
tion (RFA) techniques. These techniques include, 
for example, the well-known Rogers approximation 7 
and the Minimum State technique . 8 The RFA tech- 
niques transform frequency-domain GAFs into state- 
space (time domain) models amenable for use with 
modern control theory and optimization. This over- 
all process transforms time-domain information (CFD 
results) into the frequency-domain only to have the 
frequency-domain information transformed back into 
the time domain. Gupta et al 9 and Cowan et al 10,11 
applied a set of inputs to an unsteady CFD code and 
used the information to create an ARM A (autoregres- 
sive moving average) model that was transformed into 
state-space form. Although this technique is applied 
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entirely within the time domain, the shape of the in- 
puts applied to the CFD code requires tailoring in 
order to excite a specific frequency range, resulting in 
an iterative process. A more direct method is desired 
for this type of analysis. 

In structural dynamics, the realization of discrete- 
time state-space systems that describe the modal 
dynamics of a structure has been enabled by the 
development of algorithms such as the Eigensystem 
Realization Algorithm (ERA) 12 and the Observer 
Kalman Identification (OKID) 13 Algorithm. These 
algorithms perform state-space realizations by us- 
ing the Markov parameters (discrete-time pulse re- 
sponses) of the systems of interest. These algo- 
rithms have been combined into one package known as 
the System/Observer/Controller Identification Tool- 
box (SOCIT) 14 developed at NASA Langley Research 
Center. The present research is the first application 
of these techniques to the development of state-space 
models of linearized unsteady aerodynamic systems 
based on the modal pulse responses of a nonlinear 
(CFD) aeroelastic system. 

The method presented in this paper provides a new 
approach that bypasses the need for generating CFD- 
based, frequency-domain GAFs by directly trans- 
forming the CFD-based, time-domain modal pulse 
responses into discrete-time unsteady aerodynamic 
state-space models. These state-space models can then 
be used for computing the aeroelastic/aerodynamic 
response of the vehicle due to arbitrary motions for 
ASE analyses, including time-domain flutter analyses, 
simulations, and control-law design. The goal of this 
paper is to demonstrate the applicability of the SOCIT 
algorithm to the development of state-space models of 
unsteady aerodynamic systems using CFD-based pulse 
responses. 

This paper begins with a brief outline of the overall 
process, followed by a description of the ERA algo- 
rithm (contained within the SOCIT) and a description 
of the CFD-based pulse response technique. Details re- 
garding the computational model of the AGARD 445.6 
Aeroelastic Wing are presented including some results 
from previous research to provide the reader with an 
adequate background for the present study. Results 
are presented for state-space models generated directly 
from CFD-based pulse responses, including an assess- 
ment of the accuracy of the technique. 

Description of Methods 

Process Outline 

An outline of the process that tranforms CFD-based 
pulse responses into state-space models is as follows: 

1. Implementation of pulse (or step) response tech- 
nique into aeroelastic CFD code; 

2. Computation of pulse (or step) responses for each 
mode of an aeroelastic system using the aeroelastic 
CFD code; 


3. Pulse (or step) responses generated in Step 2 are 
the Markov parameters used by the SOCIT algorithm; 
if step responses are used, then the derivative of the 
step responses (which are pulse responses) need to be 
computed prior to implementation in the SOCIT al- 
gorithm; 

4. Evaluation/validation of the state-space models 
generated using SOCIT; this involves a comparison of 
the frequency content of the original CFD-based pulse 
responses (Step 2) with the frequency content of the 
state-space models; 

Steps 1 and 2 are described in greater detail in the 
references that address Volterra-based Reduced-Order 
Models (ROMs) such as Refs. 1-5. The basic premise 
of Volterra-based ROMs is the extraction of linear 
and nonlinear kernel functions that capture the input- 
output functional relationship between, for example, 
unsteady motion of a wing (input) and the resultant 
loads created by that motion (output). For Volterra- 
based ROMs, these kernel functions are linearized and 
nonlinear impulse response functions. When applied 
to a discrete-time system such as a CFD code, these 
impulse response functions are referred to as pulse 
response functions or, simply, pulse responses. Addi- 
tional details can be found in the stated references. For 
completeness, however, the relevant aspects of Step 2 
are discussed in the following section. 

CFD-Based Discrete Unit Pulse Response 
Technique 

Considering an aeroelastic system as the coupling of 
an unsteady aerodynamic system (CFD code, in this 
case) and a structural system (Figure 1), the present 
study focuses strictly on the unsteady aerodynamic 
model. A standard technique for computing linearized 
generalized aerodynamic forces (GAFs) for an aeroe- 
lastic system with n modes using a CFD code is the 
application of a Greens function (influence function) 
approach. Using the CFD code, each mode is individ- 
ually excited to obtain the response of all the modes to 
this excitation. This process is applied to all n modes, 
resulting in an n by n "matrix” of responses. The term 
"matrix” is in quotes to indicate that the responses 
obtained using this method are usually time-domain 
functions rather than the constants that usually pop- 
ulate a standard matrix. 

This technique is a linearization by virtue of the fact 
that, in a computational aeroelastic analysis, the in- 
put to the nonlinear flow solver is the total physical 
deformation of the wing consisting of the summed to- 
tal of all the modes of interest. By applying a separate 
excitation to each mode through the nonlinear flow 
solver, the total nonlinear aeroelastic response is being 
approximated by a linear superposition of its individ- 
ual responses. For a linear flow solver, this approach 
would be exact. In this case, because the flow solver is 
nonlinear, this approach is a linearized approximation. 
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Fig. 1 Schematic of identification of generalized 
aerodynamic forces (GAFs). 

Consistent with this assumption, this approximation 
is valid only for small input amplitudes. This is not 
necessarily a drawback as, quite often, the linearized 
dynamic aeroelastic response about a nonlinear steady 
(or static aeroelastic) condition is a reasonable repre- 
sentation of the system under investigation. 

There are three types of modal excitation inputs 
that are typically used when implementing this tech- 
nique. The first is a brute- force approach based on the 
input of sine waves of individual frequencies. The in- 
dividual modal responses to these inputs for n modes 
and r frequencies requires n times r separate code eval- 
uations. In addition, the time length required for 
each one of these evaluations can be quite large (i.e., 
computationally expensive) in order to get an ade- 
quate number of cycles for accurate post-processing, 
especially for the lower frequencies. This approach is 
clearly the least efficient. 

A second, more elegant approach, involves the use 
of an exponential (Gaussian-shaped) pulse. 6 The ex- 
ponential pulse can be shaped to excite a particular 
range of frequencies with a broad exponential pulse 
exciting primarily low frequency modes and a sharper 
exponential pulse exciting primarily higher frequency 
modes. Because an exponential pulse excites a pre- 
selected frequency range, only one code evaluation is 
required per mode. This is a significant computational 
savings compared to the brute force approach, but 
shape optimization of the exponential pulse is required 
when targeting a particular frequency range. In addi- 
tion, for low frequencies, the exponential pulse needs 
to be wide, resulting in long computational times. 

A third, recently-developed, approach consists of re- 
placing the exponential pulse input with a unit pulse 
(discrete-time equivalent of unit impulse) or step input 
consistent with the Volterra-based ROM method. 2-5 


Since the unit pulse/step input excites the entire fre- 
quency range of a system, no shape optimization is 
needed. In addition, due to the short time length 
of these responses, each code evaluation is signifi- 
cantly shorter in computational length (and cost) than 
the code evaluations for the exponential pulse input. 
Raveh et al 5 applied this technique successfully to the 
AGARD 445.6 Aeroelastic Wing using an aeroelastic 
CFD code. The shapes of the first four structural 
modes for this wing are presented in Figure 2. 



Fig. 2 AGARD 445.6 wing first four elastic mode 
shapes mapped into the CFD surface grids 

Consistent with the linearization process described 
above and in order to reduce the possibility of nu- 
merical problems with aeroelastically-deforming grids, 
small amplitudes are used with this technique. 

For the CFD computations, the flow field around 
the wing was evaluated on a C-H type grid, with 193’ 
points in the chordwise direction along the wing and 
its wake, 65 grid points in the spanwise direction, and 
41 grid points along the normal direction. The flow 
was analyzed using the EZNSS (Elastic Zonal Navier- 
Stokes Solver) Euler/Navier-Stokes code. 18 This code 
provides a choice between two implicit algorithms, 
the Beam and Warming algorithm 15 or the partially 
flux- vector splitting algorithm of Steger et al. 16 Grid 
generation and inter-grid connectivity are handled us- 
ing the Chimera approach. 17 The code was enhanced 
with an elastic capability to compute trimmed maneu- 
vers of elastic aircraft. 18 

The process of mode-by-mode excitation using var- 
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ious types of inputs, discussed previously, was per- 
formed for this wing using four elastic modes at several 
Mach numbers. The mode-by-mode excitation tech- 
nique provides the unsteady aerodynamic response in 
all four modes due to an excitation of one of the modes. 
In this fashion, the matrix of four-by-four response 
functions is developed, resulting in a total of sixteen 
response functions. 

Two sets of excitation inputs were used by Raveh 
et al: 4 5 the discrete-time unit pulse input and the 
discrete-time unit step input. Only one of these inputs 
is needed for computing the necessary responses, but 
Raveh et al 5 applied both a pulse and a step input for 
comparisons regarding numerical sensitivity of each in- 
put. The rapid convergence of the responses (to either 
a pulse or a step) resulted in computationally-efficient 
code evaluations. In addition, the mode-by-mode ex- 
citation process required a total of only four code 
evaluations (one per mode for pulse, one per mode for 
step). These sixteen pulse (or step) responses define 
the linearized Volterra-based ROMs for this configu- 
ration at this Mach number. 

Several CFD solutions due to various sinusoidal ex- 
citations of a given mode were generated, as in the 
brute force approach described above, for the purpose 
of comparison with the Volterra-based ROM approach. 
Figure 3 presents the GAFs for all four modes due to 
a 5 Hz frequency excitation of the first mode, compar- 
ing the direct CFD solutions with the results obtained 
via convolution of the modal step responses with a 5 
Hz sinusoid. As can be seen, the comparison is excel- 
lent for most of the responses with a slight discrepancy 
for the responses of the fourth mode. The full CFD 
solution, consisting of 8000 iterations required approx- 
imately 24 hours on an SGI Origin 2000 computer with 
4 CPUs. By comparison, the Volterra-based ROM 
response required about a minute of computing time 
using digital convolution. Even including the cost of 
computing the modal step (or pulse) responses, the 
computational cost savings are significant. More im- 
portantly, the same step (or pulse) functions can now 
be used to predict the response of the unsteady aero- 
dynamic system to any arbitrary input of arbitrary 
length. As additional examples, Figure 4 presents the 
resultant GAFs due to 40 and 80 Hz sinusoidal exci- 
tations of the first mode, comparing the direct CFD 
solutions with the results obtained via convolution of 
the modal step responses with the corresponding 40 
and 80 Hz sinusoids. Comparisons are excellent, ver- 
ifying the capability of the step responses to predict 
the response to arbitrary inputs. Similar results were 
obtained using the pulse responses as well. 

In order to compare with published frequency- 
domain GAFs for this configuration, Raveh et al 5 gen- 
erated the frequency-domain GAFs for the AGARD 
445.6 wing by performing convolutions of the corre- 
sponding modal pulse/step responses with various fre- 
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Fig. 3 Comparison of modal responses to a 5 Hz 
sinusoidal input for direct (CFD) and convolved 
(Volterra-based ROM) methods. 

Mode 1 Mode 2 



Non-Dimensiona! Time 

Fig. 4 Comparison of modal responses to a 40 Hz 
and an 80 Hz sinusoidal input for direct (CFD) and 
convolved (Volterra-based ROM) methods. 

quencies of interest. Another approach, as discussed 
in Ref. 3, is to perform a Fourier transform of the 
pulse responses directly, resulting in the desired full 
frequency spectrum without the need for multiple con- 
volutions. In the interest of brevity, those results 
will not be repeated here as they are discussed in 
detail in Ref. 5. However, the pulse responses com- 
puted for that study will be reviewed and discussed 
since these are the principal components used by 
the System/Observer/Controller Identification Tool- 
box (SOCIT) algorithms. 

System / Observer/Con troller Identification 
Toolbox (SOCIT) 

The primary algorithm within the SOCIT group 
of algorithms used for the present system realization 


4 OF 9 


American Institute of Aeronautics and Astronautics Paper 2001-1213 










is known as the Eigensystem Realization Algorithm 
(ERA). A brief summary of the basis of this algorithm 
follows. 

A finite dimensional, discrete- time, linear, time- 
invariant dynamical system has the state-variable 
equations 


x(k + 1) = Ax(k) + Bu(k) 

(1) 

y(k) = Cx(k) + Du{k) 

(2) 


where x is an n-dimensional state vector, u an m- 
dimensional control input, and y a p-dimensional out- 
put or measurement vector with k being the discrete 
time index. The transition matrix, A, characterizes 
the dynamics of the system. The goal of system real- 
ization is to generate constant matrices (A, B, C) such 
that the output responses of a given system due to a 
particular set of inputs is reproduced by the discrete- 
time state-space system described above. 

For the system of Eqs. (1) and (2), the time-domain 
values of the systems pulse response (discrete-time 
equivalent of impulse response) are also known as the 
Markov parameters and are defined as 

Y{k) = CA k ~ l B (3) 

with B an n by m matrix and C a p by n ma- 
trix. System realization techniques provide the con- 
stant matrices A, B, and C using Y(k). The ERA 
algorithm, 12 ” 14 similar to the Ho-Kalman procedure, 
begins by defining the generalized Hankel matrix con- 
sisting of the pulse responses (Markov parameters) for 
all input/output combinations. The algorithm then 
uses the singular value decomposition (SVD) to com- 
pute the A, B, and C matrices. Although Eq. (2) does 
not contain it, often, the direct feedthrough matrix, D, 
is required whenever the initial values of the Markov 
parameters are nonzero. 

The ERA algorithm has been used successfully for 
the identification of several experimental structural 
dynamic systems. Although the algorithm also has 
been used to extract damping and frequency infor- 
mation from CFD-generated aeroelastic transients (no 
published references), this research represents the first 
time that the ERA algorithm is applied to the devel- 
opment of unsteady aerodynamic state-space models 
using aerodynamic pulse responses (Markov parame- 
ters). The ability to generate state-space models of 
systems using pulse responses was a primary moti- 
vation for the development of aerodynamic pulse re- 
sponse functions. 3 Additional details regarding the 
ERA algorithm and its numerous applications are dis- 
cussed in the references provided. 






Fig. 5 Pulse responses due to inputs in a) Mode 
1, b) Mode2 c) Mode 3, d) Mode 4. 

Results 

CFD-Based Pulse Responses 

The results presented in this study are for a Mach 
number of 0.96. The GAF pulse responses due to pulse 
inputs in modes 1, 2, 3, and 4, respectively, are pre- 
sented in Figure 5. As can be seen, the responses 
are well-behaved and exhibit a rapid converge to zero. 
It is precisely this rapid convergence to zero that re- 
sults in an efficient computation of these functions. 
The ’’discrete” (non-smooth) nature of the responses 
is consistent with results obtained by Silva. 3 

Inspection of Figure 5 indicates that each one of the 
modes has the strongest response to its own excitation, 
a physically consistent result. In addition, those modes 
that are minimally correlated (i.e., the small response 
of mode 3 to mode 1 and the small response of mode 1 
to mode 3) can be readily identified as well. Although 
not applied in the present study, this feature of time- 
domain modal responses may prove useful in modal 
truncation and/or modal residualization studies. The 
metric used, in the present study, to determine the 
quality of an identified state-space system is the accu- 
racy to which the state-space system reproduces the 
CFD-based pulse responses presented in Figure 5. 

32nd-Order State-Space System 

The SOCIT algorithms are used to generate 
discrete-time unsteady aerodynamic state-space ma- 
trices (A, B, C, and D) using the set of modal pulse 
responses presented in Figure 5. The first system to 
be realized using this data is a 32nd order model. A 
discrete-time (z-transform) pole-zero plot of the eigen- 
values of this system is presented in Figure 6. For 
a discrete- time system, if all the poles of the system 
are contained within the complex unit circle, then the 
system is stable. This is analogous to a continuous- 
time system being stable if all of its poles are on the 
left-hand side (roots with negative real parts) of the 
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Fig. 6 Pole-zero map for 32nd order aerodynamic 
state-space system. 

complex plane. As can be seen in Figure 6, the un- 
steady aerodynamic system realized is stable, as it 
should be, with a dominant cluster of low-frequency, 
low-damped roots. 

Pulse responses for the 32nd-order state-space 
model were generated and compared with those of Fig- 
ure 5. The comparison of these responses is presented 
as Figure 7, where it can be seen that the approxi- 
mations of the sixteen GAFs are identical (to within 
plotting accuracy) to the GAFs in Figure 5. However, 
in order to thoroughly quantify the growth of error 
as the order of the model is reduced, comparisons of 
the percent error , for corresponding responses is pre- 
sented for a subset of the results. Figure 8 presents 
the percent error for the CFD and state-space modal 
pulses for the first four modes due to an input in the 
first mode. The percent error is quite small, as can 
be seen. But for the third mode (due to an input in 
the first mode), the percent error is noticeably larger. 
Although these are very acceptable values of percent 
error, the increased error for the third mode indicates 
an increased difficulty of the SOCIT algorithm to accu- 
rately capture responses that are small compared with 
the rest of the responses. This is related to the level 
of correlation between modes and can be interpreted 
as a measure of observability/controllability from one 
mode to another. This concept needs to be considered 
in the future development of these state-space models. 

Figure 9 is the magnitude and phase of the fre- 
quency response for the CFD and state-space pulses 
for the four modes due to an input in the first mode. 
As expected, the comparisons for the magnitude are 
very good for the first, second, and fourth modes with 
a slight discrepancy for the second and fourth mode 
results at very low frequencies. The error in the third 
mode manifests itself as an error within the 1 to 4 Hz 
frequency range. The phase results show good com- 
parison for the first, second, and fourth modes with 




Fig. 7 Comparison of pulse responses due to input 
in a) Mode 1, b) Mode 2, c) Mode 3, d) Mode 4 
for CFD and state-space model (32nd order). 


Percent Error Between CFD and State-Space Pulse Responses, 32nd Order 



Fig. 8 Percent error between CFD and state-space 
pulse responses, 32nd order. 

noticeable error for the third mode. 

Although a 32nd-order state-space model of the un- 
steady aerodynamic system is somewhat larger than 
might typically be desired, it is nonetheless a signifi- 
cant reduction of order when considering the alterna- 
tive: the full, three-dimensional Euler CFD computa- 
tional model and solution process. In addition, there 
is no need to perform any shape optimization of the 
pulses as these are computed only once per mode. The 
SOCIT algorithm computes state-space models very 
quickly (within seconds) on a workstation with no it- 
erative process involved unless the user is interested in 
evaluating various model orders. 

15th-Order State-Space System 

Using the SOCIT algorithm and selecting the more 
dominant components of the 32nd-order state-space 
system, a 15th-order state-space system was gener- 
ated. The z-plane pole-zero map for this system is 
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Fig. 9 Comparison of modal frequency responses 
due to input in mode 1, 32nd order. 



Fig. 10 Pole-zero map for the 15th order state- 
space model. 

presented in Figure 10 with stability of the approx- 
imated unsteady aerodynamic system preserved, as 
expected. More importantly, inspection of z-plane 
pole-zero maps are helpful in identifying the dominant 
roots of a system as the order of an approximate sys- 
tem is varied. Figure 11 presents comparisons for 
the CFD and state-space pulse responses for the 15 th- 
order model. 

The comparisons are very good but error analysis for 
the first four modes due to an input in the first mode 
indicates an increase in the error, in particular for the 
third modal response as shown in Figure 12. It is 
interesting to note that if the correlation between two 
modes is low, then the response of one of the modes 
with respect to the other has a relatively smaller am- 
plitude in comparison with the other modal responses. 
This smaller amplitude response, in turn, results in an 
increase in the error between the CFD pulse response 
and the state-space model. 




Fig. 11 Comparison of pulse responses due to in- 
put in a) Mode 1, b) Mode 2, c) Mode 3, d) Mode 
4 for CFD and state-space model (15th order). 



Fig. 12 Percent error between CFD and state- 
space pulse responses 15th order. 

Finally, Figure 13 presents the comparison of mag- 
nitude and phase for the frequency responses for the 
CFD pulse responses and the 15th-order state-space 
model. Compared with the 32nd~order model, an in- 
crease in the error for both magnitude and phase, in 
particular for the third mode, in the low frequency 
range is evident with good comparison for the re- 
mainder of the frequency range. However, there does 
appear to be an improvement in the approximation 
of the magnitude for the third mode over a range of 
low frequencies. However, the phase approximation for 
the third mode has a noticeable increase in the error. 
Additional research, currently underway, is needed to 
fully understand the impact of these modeling errors 
on aeroelastic and aeroservoelastic analyses. 

Concluding Remarks 

This study presented the direct, time-domain re- 
alization of linearized, unsteady aerodynamic state- 
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Fig. 13 Comparison of modal frequency responses 
due to input in mode 1, 15th order state-space 
model. 

space models using CFD-generated pulse responses. 
By performing this realization entirely in the time do- 
main, the traditional and inefficient approach of trans- 
forming time-domain GAFs into the frequency domain 
only to convert these back into the time domain via ra- 
tional function approximations was avoided. 

The results presented were for the AGARD 445.6 
Aeroelastic Wing using pulse responses that were com- 
puted for this wing using the EZNSS Euler/Navier- 
Stokes flow solver with aeroelastic capability. The 
pulses were computed for a transonic Mach number 
of 0.96 and a mode-by-mode excitation process re- 
sulted in a linearized set of GAF influence functions. 
The System/Observer/Controller Identification Tool- 
box (SOCIT) was then used to generate four-input, 
four-output discrete-time state-space models of vari- 
ous orders that accurately capture the dynamics of 
the CFD-generated unsteady aerodynamic responses. 
A 32nd-order state-space model was shown to capture 
the frequency content of the pulse responses with neg- 
ligible errors. Further model-order reduction resulted 
in a 15th-order state-space model with slight increases 
in the error for low frequencies but with excellent cor- 
relation for the higher frequency range. 

The results provided some insights into the appli- 
cability and efficiency of the technique. It was no- 
ticed that if a given mode does not strongly excite 
another mode, then that portion of the identification 
process results in larger errors. This appears to be re- 
lated to an observable/controllable aspect of the modal 
system but additional research is needed to fully un- 
derstand this effect. Future research will address the 
inclusion of these state-space models into an aeroelas- 
tic/aeroservoelastic analysis to better understand the 
impact of these mathematical approximations. Com- 
parisons of these systems with the more traditional 
rational function approximations will be performed as 


well. Finally, the results presented provide a first step 
towards the development of bilinear state-space ma- 
trices that will incorporate the effects of the nonlinear 
(higher-order) kernels. 
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